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Abstract 

When it is polarised, a cell develops an asymmetric distribution of 
specific molecular markers, cytoskeleton and cell membrane shape. Po- 
larisation can occur spontaneously or be triggered by external signals, 
like gradients of signalling molecules... In this work, we use the models 
of cell polarisation introduced in [9] and [4] and we set a numerical anal- 
ysis for these models. They are based on nonlinear convection-diffusion 
equations and the nonlinearity in the transport term expresses the pos- 
itive loop between the level of protein concentration localised in a small 
area of the cell membrane and the number of new proteins that will be 
convected to the same area. We perform numerical simulations and we 
illustrate that these models are rich enough to describe the apparition 
of a polarisome. 

Keywords: Cell polarisation, global existence, blow-up, numerical 
simulations, Keller-Segel system. 

1 Introduction 

Cell polarisation is a symmetry-breaking event that occurs in cell division, 
mating or morphogenesis. Molecular markers play a central role in estab- 
lishing this phenomenon. Indeed, there are two different behaviours: a 
non-polarised cell has its markers radially homogeneously distributed while 
markers are located in a small area of the cell membrane for a polarised 
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cell. Yeast cells are dynamically polarised in response to the extracellular 
gradients of signals (chemokynes). However, it has been observed in [TJ] 
that polarisation can occur spontaneously without any external asymmetric 
stimulus. 

During the past decade, many models describing cell polarisation have 
been developed. The majority of these models are based on reaction-diffusion 
systems where polarisation appears as a type of Turing instability [10], [12] . 
[TT] , or due to stochastic fluctuations [3] , other models include cytoskeleton 
proteins as a regulatory factor [7J, [14] . Many biological studies have shown 
that the cytoskeleton plays an important role in polarisation. It has been 
suggested that there is a positive feedback on molecular markers density. 
Indeed, disruption of transport along the cytoskeleton greatly reduces the 
stability of polar cap [14] . The cell cytoskeleton is a network of long semi- 
flexible filaments made up of protein subunits [13] . These filaments (mainly 
actin or microtubules) act as roads along which motor proteins are able to 
perform a biased ballistic motion and carry various molecules. Molecular 
markers play a key role in the formation of these filaments. 

Following [9], [5] and [I], in this work we study models that describe 
the dynamics of cell polarisation. In these models, molecular markers, such 
as proteins, diffuse in the cytoplasm and are actively transported along the 
cytoskeleton. The resulting motion is a biased diffusion regulated by the 
markers themselves. Using numerical simulations and mathematical heuris- 
tics, we observe that the coupling on the velocity field achieves an inhomo- 
geneous distribution of molecular markers without any external asymmetric 
field. Such an inhomogeneous distribution is only due to interaction between 
molecular markers. 

Throughout this paper, the density of molecular markers (resp. advec- 
tion field) is denoted by p(t, x) (resp. u(t, x)). The advection is obtained 
through a coupling with the membrane concentration of markers. The cell 
is figured by the domain Q C R n with n = 1,2 and a part of the boundary of 
the domain will be the active membrane denoted by T. The time evolution 
of the molecular markers satisfies the following advection-diffusion equation, 
see [9] and [1]: 



There is no creation nor degradation of molecular markers in the cell, so the 
quantity of molecular markers remains constant in time: 



This condition is ensured by a zero flux boundary condition on the boundary. 
A first simplified step is to assume that the cell is essentially bidimensional 



%>(i,x) =DAp(f,x)-xV.((3(t,x)u(t,x)), t>0, xe!l, 
p(0,x) = p (x). 
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and to neglect curvature effects. The membrane boundary is then a ID line 
along the y-axis and the cytoplasm is parametrized by x = {x, y) £ M + x R. 

The plan of this work is the following. First, we present the models, that 
are based on different expressions for u, and we recall the main mathemati- 
cal results of the simplified model in ID for = (0, oo) and r = {x = 0}, see 
[5], [3] for more details. Then we study a more realistic model, that includes 
dynamical exchange of markers on the boundary. This model was intro- 
duced in [9] and studied in jlj. We provide a methodology for parameter 
estimation and qualitative description of cell polarisation by using mathe- 
matical heuristics. Then, we perform a numerical analysis of this model. We 
introduce the numerical part by the one dimensional case. Finally, we give 
tools to study the numerical implementation of the model on an annulus 
domain. 

2 Presentation of the models and mathematical 
results 

2.1 One dimensional case 

In this section, we study the one dimensional case on the half line for f2 = 
(0, oo). The membrane is then the point T = {x = 0}. For the first model, 
the advection field towards the membrane is equal to the density of molecular 
markers on the boundary p(t, 0). Then we improve this model by considering 
that only the trapped molecular markers on the membrane contribute to the 
advection field. 

2.1.1 Simplified model set on the half line 

In [3] a first mathematical studies has been done on this model. We define 
an advection field u(t, x) for 

u(t,x) = -p(t,0), 

in such a case reads as (with D = 1 and x = 1): 

d t p(t,x) = d xx p(t,x) + p{t,0) d x p(t,x), t>0, x > 0, (3) 

with the following zero flux condition on the boundary T = {x = 0}, that 
ensures the mass conversation ([2]), 

d xP (t, 0) + p(t, 0) 2 = 0. (4) 

Interestingly enough in [1] , it has been proved that solutions of ([3| blow-up 
in finite time if their masses are above a certain critical mass, M > 1, and 
exist globally in time if M < 1. Let us first recall the definition of weak 
solutions of ([3]). 



3 



Definition 2.1 We say that p(t,x) is a weak solution of ^ on (0, T) if it 

satisfies: 

p G L°°(0, T; L]j_(M + )) , d s p € L\(0,T) xR+) , 

and p(t,x) is a solution of ^ in the sense of distributions in P'(R+). 

Let us now recall the main results for weak solutions of ([3]). 

Theorem 2.2 (Global existence: M < 1) Assume that the initial data 
po satisfies both po G L l {(\ + x) dx) and f x>Q po(x)(log po(x)) + dx < +oo. 
Assume in addition that M < 1, then there exists a global weak solution of 
equation 

Theorem 2.3 (Blow-up: M > 1) Assume M > I. Any weak solution of 
equation ^ with non-increasing initial data po blows-up in finite time. 

Remark 1 It would tempted to interpret blow-up of solutions of the one 
dimensional model as cell polarisation but concentration of markers on the 
boundary doesn't automatically mean polarisation. Indeed, consider a radi- 
ally symmetric 2D cell case then equation ([T]) reduces to the one dimensional 
one. Above a threshold on the total mass, the convection wins and mark- 
ers concentrate on the boundary. In some situations, these markers may be 
homogeneously distributed on the boundary and in such a case there is no 
symmetry breaking. 



2.1.2 The model with dynamical exchange of markers at the 
boundary 

Such a direct activation of transport on the boundary seems to be unrealistic. 
Indeed possible occurrence of blow-up in finite time suggests this claim. We 
improve the previous model by distinguishing between cytoplasmic content 
p(t, x) and the concentration of trapped molecules on the boundary, that will 
be denoted by p(t). The dynamical exchange of markers at the boundary 
is done with an attachment rate k on and a detachment rate k ff, hence the 
time evolution of p(t) is 

4/4*) = kon p(t, 0) - k ff p{t). (5) 

The advection field u(i, x) in ([!]) is now defined by 

u(t,x) = -p(t), 
hence ([!]) (with D = 1 and x = 1) reads as: 

d t p(t,x) = d xx p(t,x) + p(t)d x p(t,x), t>0, x>0, (6) 
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with a modified boundary condition 



d xP (t,0) + p(t,0)v(t) = ~»(t). (7) 
This ensures the following mass conservation shared among p(t, x) and fi(t): 

M = I po(x)dx + fio = / p(t,x)dx + fi(t). (8) 
Jr+ ' Jr + 

With equation the self-activation of transport by p(t, 0) is then delayed 
in time. Since the transport speed is bounded p,(t) < M, the solution of the 
model with dynamical exchange on the boundary exists globally in time. 
More precisely it is possible, see @], to prove that it converges towards a 
non trivial stationary state. 

Theorem 2.4 (Global existence: dynamical exchange case) Assume 

that the initial data po satisfies both po G L 1 ((l+x) dx) and f x>Q po(x)(\og po(x)) + dx < 

+oo. Assume the mass is super- critical M > 1. The partial mass m(t) = 

Ix>o Pty> x ) ^ x conver 9 es t° 1 an d the density p(t, x) strongly converges in 

L 1 towards the exponential profile (M — 1) exp(— (M — i)x). 



2.2 Two dimensional case 

In this section, we study the two dimensional case on n C M 2 . As in the 
one dimensional case, the advection field towards the membrane depends on 
the density of molecular markers on the boundary, two different situations 
(actin and microtubules) will be described, then we improve this model by 
considering an exchange of markers at the membrane and only the trapped 
molecular markers contribute to the advection field. 



2.2.1 Simplified model set on the half plan 

We study the model on the half plan for x = (x, y) £ = R + x R. The 
membrane is then the line r = {x = 0} x R. We have the following boundary 
condition for p(t, x) at point x G T 

(D Vp(t, x) - X p(t, x) u(t, x)).n x = 0, (9) 

where n x is the outward normal to T. This ensures the following mass 
conservation: 

M= p (x)dx= / p(t,x)dx. (10) 
Jn isj 

First we consider the transversal case, the field u is normal to the boundary 

u(t,x) = -S(y)p{t,0,y)e x . (11) 
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The microtubules of the cytoskeleton are normally oriented to the cell mem- 
brane and their growth depends on the density of molecular markers on 
the boundary. In the potential case, we consider the following advection 
field deriving from a harmonic potential modelling the transport by actin 
filaments: 



u(t, x) = Vc(i,x), where 



-Ac(t,x) = 0, 



Vc(t, x).n x = S(x)p(t,x 
This advection field orientation is due to the actin networks. 



if x G O, 

if x e r. 



(12) 




Actin filaments are attached on the membrane and randomly distributed, 
there orientations are mixed up. We also add the external pheromone con- 
centration at x G r which acts by the mating-pheromone MAPK cascade 
on the actin transport. We have global existence and blow-up theorems for 
the simplified model with x = (x, y) G Vt = M + x M. For clarity, we recall 
this result, see [3] for more details. 

Theorem 2.5 (Global existence in dimension 2) Assume that the ad- 
vection field satisfies the two following conditions: V • u > and u(t, 0, y) ■ 
e x = —p(t,0,y). Assume that the initial data po satisfies both po 6 L ((1 + 
|x| 2 )dx) and 1 1 Po 1 1 z^ 2 * s smaller than some constant c. Then there exists a 
global weak solution to equations Q-Q. 

Theorem 2.6 (Blow-up in dimension 2) Assume that p(t,x) is a strong 
solution to ([I]) which verifies: 

• d x p(t, x) < for all x £ Q and t > when the advective field is given 



by (11). 



d x p(t, x) < for all x G SI and t > the matrix A(t, x) = x (g) 
d x d y log p(t, x) satisfies A T + A > (in the matrix sense) when the 



advective field is given by (12) 



Assume in addition that the second momentum is initially small enough: 
there exists a constant C such that J xg Q |x| 2 /9o(x)dx < CM 3 . Then the 
maximal time of existence of the solution is finite. 
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2.2.2 The model with dynamical exchange of markers at the 
boundary 

Let 1] C t 2 be the cytoplasm domain, as in the one dimensional case ^ we 
consider dynamical exchange of markers at the boundary, so for x £ T we 
have the evolution in time of p(t, x) 

d t n(t, x) = k on p(t, x) - k ff p(t, x), (13) 

with a modified boundary condition for p(t, x) at point x G T 

(DVp{t,x) - xp(i,x)u(t,x)).n x = -d^(t,x), (14) 

where n x is the outward normal to T. This ensures the following mass 
conservation sharing by p(t, x) and p(t,x): 

M = I po(x)<ix+ / po(x.)dx = / p(i,x)cbc+ / /i(t,x)cbc. (15) 
Jn Jr Jn Jr 

As before, we consider the following advection field (cytoskeleton): 

r+ ^ v f+ \ ^ J-Ac(t,x) = 0, if x eft, 

u(t,x) = Vc(i,x), where < o/ w + ^ - f c r (16) 

I Vc(i,x).n x = o(xJ/i(r, xj, it x G 1 . 

We can also consider the transversal case. For example for ft = IR + x M, we 
take similarly as before 

u(t,x) = -S(y)fi{t,y)e x . (17) 

For the model with exchange on the boundary, blow-up or global existence 
have not been proved yet. In this work, we make a first step in this direction 
by using a mathematical heuristic and numerical simulations. There is also 
one open question: does advection field (transversal or potential) create a 
break of symmetry ? 



3 Heuristic 

The mathematical analysis performed in [1] has demonstrated that a class 
of models exhibit pattern formation (either blow-up or convergence towards 
a non homogeneous steady state) under some conditions. However the main 
question still remains unanswered: do these models describe cell polarisation 
or not? Thus in order to provide a first answer to this question, in the next 
section, we will perform numerical simulations. Our aim will be to see if, 
under some conditions, the model leads to a concentration of markers, not 
only on the boundary, but on a small region of the boundary. 

Let us now describe the mathematical heuristic which is done on the half 
plan ft = M + x IR. We give formal arguments to motivate the differences 
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arising in the dynamics leading by the two possible drifts uy (transversal 



case) and up (potential case). Let up = Vc be the solution of (12) we have 
(see 0) 

c(x, y) = - 1 / \og^{y-y'Y + x^)S{y') P {t, 0, y')dy' . 



We notice that the two possible drifts up and up share common features: 
they are both divergence free and their normal components at the boundary 
coincide. 

up • e x = up • e x = -S(y)p(t,0,y) . 

On the other hand, a key difference holds when looking at the tangential 
component at the boundary: 

up • e y = , 

u P -e y = -itH(S(-)p(t,0,-)), 
where T~L denotes the one-dimensional Hilbert transform: 



W)(y) = -p-v. / —f(x)dx. 



We expect the solution to concentrate on the boundary in the super-critical 
case for both choices of u, numerical simulations suggest it (see section [473] ) . 
Postulating the ansatz p(t,x,y) = v{t,y)8{x = 0), we can formally write 
the dynamics of v (t, y) for the two cases. Integrating the main equation ([!]) 
with respect to x with zero flux condition on V = {x = 0}, we obtain: 

d t I p(t,x,y)dx = Ddyy [I p(t,x,y)dx)-xd y (/ p(t,x,y)(u(t,x,y) ■ e y )dx 
Jr + \Jr + J \Jr + 

• up • e y = 0: diffusion equation for v: 

d t u(t,y) = Dd yy u(t,y) . 

• up • e y = —n7i(S(-)p(t, 0, •)) and assuming S constant on R, it reads 
as: 

d t u{t,y) = Dd yy v{t 1 y)+xSdy{v{t,y)U{v){y)) . 

Transversal dynamics are very different: in transversal case boundary dif- 
fusion dominates while in potential case Hilbert transform has a critical 
singularity to offset the diffusion on this equation. This latter equation ex- 
hibits blow-up if f R v(t,y) dy = M is above the critical mass see [6J 
done in a peculiar variant of Keller-Segel equation for more informations. 
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Remark 2 This is a first step to observe a critical mass phenomenon and 
this may lead to blow-up if the mass is large enough. In this way, we also 
define an order of magnitude for some parameters. It is to be noticed that 
this latter criterion is valid for an infinite domain, namely y £ R. In the 
case of a cell, the domain is finite and the existence of such a dichotomy 
has not been proved yet. In order to obtain more information on the critical 
value distinguishing the polarised case and the stable case we will perform 
numerical simulations. 

This heuristic shows us that potential case is able to break symmetry 
more readily than the transversal case. This difference between the two 
models is also discussed in the following section. 



4 Numerical analysis 

We introduce this numerical section by the discretization of the convection- 
diffusion model set on a ID periodic domain. This first step allows us 
introducing the discretization of this model on a 2D domain. We model the 
cell as an annulus, molecular markers cannot pass inside the nucleus. In this 
section, for simplicity we fix all the parameters values to 1 except M. 



4.1 First step with the one dimensional case 

Let u(t,9) be a periodic function on the periodic domain £1 = M/27rZ. We 
consider the following advection-diffusion equation 

d t p(t,6) = d e {d e p{t,6) - p{t,9)u(t,9)), t > 0, 9£R/2ttZ. (18) 

Let t n = n At be the time discretization and {9 k = k A9,k £ {1,..., Ng } } be 
the space discretization of the periodic interval R/27rZ. Since the equations 
of the model are written in a conservative form, the natural framework to 
be used for the spatial discretization is the finite volume framework. We 
hence introduce the control volume defined for k £ {1, ...,Ng} 

V k = (9 k _i,9 k+i ). (19) 



2 1 2 



Let p k (resp. u? i ) be the approximated value of the exact solution p(t n , 9k) 
(resp. u(t n ,9 k+ i)), the classical upwind scheme for (18) reads as 

D n + l _ Q n F k , i - T k _ i 

Pj ^T K = \ e - > k£{l,...,N 9 }, (20) 
where the numerical flux T k+ i and i r fe _i are defined by 
n+l _ n+1 

X- , _ yk _ A up ( n+l n+l n+U 



2 



2 



n+l _ n+l 
_ P k Pk-l _ A u p( n+l n+l n+U 
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with the advection numerical flux is given by 



A up (u,x„,x+) 



UX- 



ux< 



if u > 0, 
if u < 0. 



(21) 



The periodic flux condition on boundary reads as 



2 



n+1 



AO 



and we recall that u is periodic so we set the value u™ = -u™ , 1 . The diffusion 
and advection terms are both treated implicitly, the scheme is then uncon- 

rp 

ditionally stable. We define the column vector p n = (p™ p^ ... PAr e ) • 
As usual, see e.g. [2], the discrete heat matrix A G MN g (M) with periodic 
flux condition on the boundary is defined as 



/2 



A 



1 



V- 1 



(22) 



/ 



Periodic flux condition adds the top right term and the bottom left term. 



Next, in order to use A up defined by equation (21), we define (u) + = 
max(it, 0) and = min(u,0) so A up (u, p k , Pk+i) = {u) + Pk + (u)~ Pk+i- 
The discrete advection matrix B n £ Mn 6 (M) with periodic flux condition 
on the boundary is then defined as in [2] 



B' 



u n 3 



U N e + \ 



\ 



(23) 



u n x 



J 2 



Ne-h) ) 
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At each time step we have 

n+l _ n i I 

At ~ A9* Ap A9 B 9 ■ 

We use a standard numerical method to invert the matrix A + A9 B n+l + 
-£rlNe- Finally, at each time step we resolve 

P ^ = (A + A6B^ + ^-I N y A ° 2 



At JV V At 
4.2 Two dimensional case: polar case 

Let Q, C M 2 be the cytoplasm domain, the domain is obviously an annulus 
where the molecular markers cannot enter in the nucleus of the cell. Let us 
first recall the model on the annulus ft = B(0, R max ) \ B(0,R m i n ) with n x 
the unit normal vector to f2 at point x G dQ (we note C(0, R) the circle of 
center (0,0) and radius R): 

d t p(t, x) = V. (Vp(t, x) - p(t, x) u(i, x)) , in ft, (24) 
(Vp(t, x) - p(t, x) u(t, x)).n x = 0, on C(0, i? maa: ), (25) 
(Vp(t, x) - p{t, x) u(t, x)).n x = 0, on C(0, # min ). (26) 

Since the domain Q is assumed to be an annulus, it is appropriate to intro- 
duce polar coordinates r and 9. Let x = (r cos(#), r sin(#)) G we have 
the following equations on ^p(t,r,9) = p(t,x) with (r,9) G [i?mm, fimai] x 
M/27TZ: 

$p(t, r, 0) = d r (rd r ( P^ r '°A - p( t , r, 9)u r (t, r, 9) 



+ »« ( ^(d e p(t,r,9) - p(t,r,9)u d (t,r,9)) ) . in O. (27) 



I) ri), ( P^llS. )- f) (l. r.(~))u r (l.r.O). on ('( 0. . /?,„,,,.).( 28 ) 

= r9 r ( ^ t,r ' g) ) -p(t,r,9)u r (t,r,9), on C(0, i? mm ).(29) 



We recall the potential case u(i,x) = Vc(i,x). Laplace equation on c with 
non appropriate Neumann conditions on a bounded domain is ill-posed, see 
[2] e.g. In order to handle this problem, we add a degradation term: 

-Ac(x) + ac(x) = 0, in (30) 
Vc(x).n x = p(i,x), on (7(0,^^), (31) 
Vc(x).n x = 0, on C(0, i? min ). (32) 
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We also give the previous equations in polar coordinates by -c(r, 9) = c(x) 
With (r,0) G [i?™n,^max] x M/27rZ: 

(r0 r (^)) - ^%c(r,e) + ac(r,9) = 0, in ft, (33) 
dr ^) = on C (0,R max ), (34) 

S r (^) =0, on<?(0,iW)- (35) 

If we consider dynamical exchange of markers at the active boundary, for 
x € T = C(0, Rmax) we have the evolution in time of p(t,x). 

d t p(t,x) = p(t,x) - yu(i,x), on C(0, i? max ). (36) 

We replace then p8| by 



fyz(*>0) = rd r ( P(t,r,g) ) -p(i,r,0)u r (i,r,0), on C(0, R max ){37) 



and (34) by 



5r ( ) = e), on C(0, iW). (38) 



r 

The transversal case is then 

u(t,r,6) = e r , 

"mm 

and we can also write the transversal case for the dynamical exchange model 

u(i,r,0) = p(t,9)e r , 

Let t n = n At be the time discretization and {rj = R m %n + j Ar, j € 
{1, iVj.}} (resp. {6>& = kA9,k G {1, iVp}}) be the space discretization 
of the bounded interval [i? m m,i? ma i] (resp. periodic interval M/27rZ). We 
introduce the control volume IV^ u C M 2 



W (ilfc) = (r._i,r. +| )x(^_i,^ +| ) 



Let P,™ (resp. ^) be the approximated value of the exact solution 
p(t n ,rj,9k) (resp. p(t n ,9k))- Let cq-.^) be the approximated value of the 
exact solution c(rj,9k)- 

4.2.1 Equation on p 

In the dynamical exchange model, we can resolve at each time step the 



discretization of equation (36) for k G {1, N y } 
p n k +1 =pt + At(pl-pl). 
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4.2.2 Equation on c 

For simplicity, we call T the numerical flux as in the ID case, we can write 



the following scheme for equation (33) for (j, k) £ {1, N r } x {1, Ng} 



'T, 



UH> k ) 



Ar 



+ 



AO 



+ ac, 



In order to use finite volume we define the numerical flux by 



C (j + 1,*0 _ C (j,fc) 
r j + l r j 

Ar 

1 cq.fc+1) ~ C(j-,fc) 



r „• i 



c (j,fc) _ c o-i,fc) 



(i 2' fc ) ' ■> i _\.r 



(?.*-§) 



1 d U,k) ~ %,k-l) 



AO 



The zero flux boundary condition (35) impose that Tn 



and the 



boundary condition (34) F^ Nr+ i k ^ = r Nr+ i /i£ for k E {1, Ng}. Similarly, 



the periodic conditions impose for j € {1, N r } 



1 5 faM) ~ d U,N e ) 
r? AO 



We define the column vector C by C(k + (j — l)Ng) = C(j.k) with (j,k) £ 
{l,...,N r }x{l,...,Ng}; 

C = (5(1,1) • • • C(1,JV,) C(2,l) • • • d (2,N e ) ■ ■ ■ C(N r ,N e )) T 

For Ar = AO the rigidity matrix A is defined by 

/ \ 



A 



r . i r. i+r.,1 



7d 



r j+ i 



Id 



+ 



\.A 



N r -1 



A 



V 



-A 



(39) 



The flux boundary condition C(0, R ma x) imposes this right hand side column 
vector of length N r Ng: 

( o \ 



N r 







or in the exchange case 7£™ = r N , i 



/ \ 
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We use a standard numerical method to invert the symmetric positive defi- 
nite matrix -^piA + alN r N g and then resolve at each time step 

i V 1 i 

A + aI NrN A —K n c . 



^Ar 2 r 9 ) Ar 

4.2.3 Equation on p 

For simplicity, we call J- the numerical flux as in the ID case, we can write 



the following scheme for equation (27): for (J, k) 6 {1, ...,N r } x {1, ...,Ng} 



jjS) ~ %fc) = J7 u+b k ) • F UzM ^j+j) ^fc-D 

At Ar A0 

We define the numerical flux by 

pn+l pn + l 

X* - t- , r J+ 1 r o _ Aup I „.n+l pn+l pn+l \ 

^u+hk) - r j+\ Ar A X^u+hky ^ky-^{j+i,k)J ' 

pn+l pn+l 

x* 1 — rp -. r J- 1 A U P i, n+1 pn+l pn+l \ 

■^0-3.*) ~ r i-| Ar ^(j-ifej'-ny-i.*)'^)^ ' 

/ pn+l _ pn+l \ 
. 1 0>+l) 0» 4 up / ,,n+l pn+l pn+l \ \ 

/pn+l _ pn+l , \ 

77 _ _ | Q' fc ) (■?>-!) _ /i up / „n+l pn+l pn+l 

In the transversal case ( 11 ), we take u r = ^( R ™ a3; ' e ) anc i u# = we define at 
time t n 



u?. , i , s = ( jVr ' fc ) , or j n the exchange case iff. , i = — ul! - , 

0+2 > fc ) rjv r 0+2. fc ) fe 

it?- i ,\ = Nr or in the exchange case u?. i = — it?, 

0-2' fc ) r N 0-2. fc ) ^ fe 



w 0,*+§) = °> u u*-h) = °- 



In the potential case (12), we have u r = d r (^) and ug = ^(%c we define at 
time t n 

ti"_. i , s = — t — , u, . i 



0+3- fc ) Ar ' 0-a.*) Ar 

1 C(j,fc+1) ~ C(j,fc) n 1 5 0,fc) ~ g 0,fc~i) 

0>+§) rj A# ' n 0.fc-|) n AO 
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The zero flux boundary condition (29) impose that J-^i ^ =0. In the sim- 



plified model, the boundary condition (28) T, 



"(iV r +ik) 



for k G {l,...,No} 



and in the model with exchange, we have F^ Nr+ i k ^ 



n + l_ 



■ " k At H for k G 



{1, ...,Nq}. Similarly, the periodic conditions impose for j G {1, N r } 

jn+I _ pu+1 . 

_ Aup I ,.n+l pn+1 pn+1 



7; 



/ pn+i _ pn-\ 

1 \ a,^) 



We define the column vector V n by + (j — l)Ng) = PV- ^ with (j, k) G 
{l,...,iV r } x {l,...,iV e }: 



pn pn pn pn pn 

^(1,1) ' • • (1,N$) ^(2,1) • • • ^(2,7V fl ) • • ' ^(N r ,N e ) 



We define the following diagonal matrices for j G {1, N r }, U. i G M/v e (l 
and 177 i G Mjv fl (R): 

■7+2 



V 

(■ 



( u 0-+i, fe -i))" 



(u ?;+§,^ 



(3 + 5.*+l) 



( u tf+§ )fe +i))" 
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Thus we can define: 



(Ut U. 

2 



u + 



u: 



3+2 J+5 



u 



N r -\ 
J 



( o 

ut 

2 



3~2 3~2 



u 



N r 



+ 



' N r -1 



B n 



V 



r N r 



In the simplified model, we have at each time step 
-pn+i _ -pn i i 

- — ^- = --^AV n+1 - ^B n+1 V n+1 . 

At Ar 2 Ar 

In the exchange model, the flux boundary condition on (7(0, R ma x) imposes 
this right hand side column vector of length N r Ng: 



( 



'^exchange 



\ 






// n + 1 — ll n 



We have at each time step 



At 



Ar 



AP 



n+l 



Ar 



E 



— — v n 

^ '^exchange' 



We use a standard numerical method to invert the matrix A + Ar B n+1 + 
lN r N B and then resolve at each time step 



V n+1 = ( A + ArB 



n+l 



Ar 2 
~At 



I N r N e 



Ar 2 
At 



exchange 



4.3 Graphics 

We use the numerical analysis done in this article to implement it with 
Matlab. Simulations have first been done in the transversal case u-p- The 
following behaviour was obtained: 
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Figure 1: Numerical Simulations on Q = [0.2, 2.5] x M/27rZ and all parame- 
ters equal to 1 with random initial conditions. For M = 10 greater enough, 
no symmetry breaking appears, molecular markers are uniformly distributed 
on the membrane. 



The dichotomy planned on M by the heuristic holds true. Indeed we have 
done simulations for small M. 




Figure 2: Numerical Simulations on = [0.2,2.5] x M./2ttZ and all param- 
eters equal to 1 with random initial conditions. Fo. For M = 0.01 small, 
steady state is isotropic. 



Then, simulations have been done in the potential case up: 
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Figure 3: Numerical Simulations on £1 = [0.2, 2.5] x M/27rZ and all parame- 
ters equal to 1 with random initial conditions. For M = 10 greater enough, 
symmetry breaking appears. Molecular markers are concentrated on one 
point of the membrane in finite time. 



The heuristic done in section [3] allowed us thinking that the potential case 
could break symmetry more readily than the transversal case. We have also 
done numerical simulation in potential case for small M and we found that 
the numerical behaviours are similar in the two possible drifts, see Fig [2] and 
II 
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Figure 4: Numerical Simulations on = [0.2,2.5] x M./2ttZ and all param- 
eters equal to 1 with random initial conditions. Fo. For M = 0.01 small, 
steady state is isotropic. 



We have assumed that the cell occupies a circle of radius R > 0. Furthermore 
for simplicity, we consider a bounded-periodic domain 0, = [0, R] x M./2'kR'L. 
With the numerical analysis done in [Tj for the potential case, we see that 
the behaviours are similar with the annulus case developed in this article. 

5 Conclusion 

Polar description of the cell has been described in this article, this improve- 
ment fitted the real cell shape and was a first step before establishing a 
model for membrane deformation. In this work we have provided a first an- 
swer to the following question: do the nonlinear convection-diffusion models 
given in [9j and |4J describe cell polarisation or not? To do so we have used 
both a mathematical heuristic and numerical simulations, which have en- 
sured us that solutions develop symmetry breaking over a critical value M*. 
This has given us a first justification of the mathematical heuristic. On this 
point, the numerical behaviours are close to cell behaviours during biological 
experiences. In order to fit biological measurements, the choice of parame- 
ters is essential and we refer to biological literature. Several measurements 
on polarisation time and localisation of polar cap have been made, we will 
describe them in a further work. 

Acknowledgement: this research has been supported by ANR program JCJC 
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for helpful discussions. 
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